Moving from Mayavi to Matplotlib

A long standing dependency of MotionClouds is MayaVi. While powerful, it is tedious to compile and may discourage new users. We are trying here to show some attempts to do the same with matplotlib.

In [2]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

Testing 3D plots in pylab

In [3]:
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
cset = ax.contour(X, Y, Z, cmap=cm.coolwarm)
ax.clabel(cset, fontsize=9, inline=1)
  Hide output
In [4]:
X, Y, Z
  Hide output
Out[4]:
(array([[-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        ..., 
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5]]),
 array([[-30. , -30. , -30. , ..., -30. , -30. , -30. ],
        [-29.5, -29.5, -29.5, ..., -29.5, -29.5, -29.5],
        [-29. , -29. , -29. , ..., -29. , -29. , -29. ],
        ..., 
        [ 28.5,  28.5,  28.5, ...,  28.5,  28.5,  28.5],
        [ 29. ,  29. ,  29. , ...,  29. ,  29. ,  29. ],
        [ 29.5,  29.5,  29.5, ...,  29.5,  29.5,  29.5]]),
 array([[-0.01 , -0.011, -0.013, ..., -0.015, -0.013, -0.011],
        [-0.011, -0.013, -0.015, ..., -0.018, -0.015, -0.013],
        [-0.013, -0.015, -0.018, ..., -0.02 , -0.018, -0.015],
        ..., 
        [-0.012, -0.014, -0.017, ...,  0.029,  0.03 ,  0.031],
        [-0.011, -0.013, -0.015, ...,  0.016,  0.017,  0.018],
        [-0.01 , -0.012, -0.014, ...,  0.007,  0.008,  0.009]]))
In [6]:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np

n_angles = 36
n_radii = 8

# An array of radii
# Does not include radius r=0, this is to eliminate duplicate points
radii = np.linspace(0.125, 1.0, n_radii)

# An array of angles
angles = np.linspace(0, 2*np.pi, n_angles, endpoint=False)

# Repeat all angles for each radius
angles = np.repeat(angles[...,np.newaxis], n_radii, axis=1)

# Convert polar (radii, angles) coords to cartesian (x, y) coords
# (0, 0) is added here. There are no duplicate points in the (x, y) plane
x = np.append(0, (radii*np.cos(angles)).flatten())
y = np.append(0, (radii*np.sin(angles)).flatten())

# Pringle surface
z = np.sin(-x*y)

fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot_trisurf(x, y, z, cmap=cm.jet, linewidth=0.2, alpha=.6)
  Hide output
Out[6]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x10a16d590>
In [10]:
"""
Contour plots of unstructured triangular grids.
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as tri
import numpy as np
import math

# First create the x and y coordinates of the points.
n_angles = 48
n_radii = 8
min_radius = 0.25
radii = np.linspace(min_radius, 0.95, n_radii)

angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
angles = np.repeat(angles[...,np.newaxis], n_radii, axis=1)
angles[:,1::2] += math.pi/n_angles

x = (radii*np.cos(angles)).flatten()
y = (radii*np.sin(angles)).flatten()
z = (np.cos(radii)*np.cos(angles*3.0)).flatten()

# Create a custom triangulation
triang = tri.Triangulation(x, y)

# Mask off unwanted triangles.
#xmid = x[triang.triangles].mean(axis=1)
#ymid = y[triang.triangles].mean(axis=1)
#mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
#triang.set_mask(mask)

plt.figure()
plt.gca(projection='3d')
plt.tricontour(triang, z)
  Hide output
Out[10]:
<matplotlib.tri.tricontour.TriContourSet instance at 0x10bdd23f8>
In [12]:
from __future__ import print_function
"""
A very simple 'animation' of a 3D plot
"""
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import time

def generate(X, Y, phi):
    R = 1 - np.sqrt(X**2 + Y**2)
    return np.cos(2 * np.pi * X + phi) * R

plt.ion()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

xs = np.linspace(-1, 1, 50)
ys = np.linspace(-1, 1, 50)
X, Y = np.meshgrid(xs, ys)
Z = generate(X, Y, 0.0)

wframe = None
tstart = time.time()
for phi in np.linspace(0, 360 / 2 / np.pi, 100):

    oldcol = wframe

    Z = generate(X, Y, phi)
    wframe = ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2)

    # Remove old line collection before drawing
    if oldcol is not None:
        ax.collections.remove(oldcol)

    plt.pause(.001)

print ('FPS: %f' % (100 / (time.time() - tstart)))
  Hide output
FPS: 175.698996

However, there is nothing really easy yet to produce a contour plot of a scalar field, as mayavi does very well.

Testing the utility functions of Motion Clouds

Motion Clouds utilities

Here, we test some of the utilities that are delivered with the MotionClouds package.

In [1]:
import MotionClouds as mc
  Hide output

progressbar

Wrapper around pyprind.

In [2]:
### PERCENTAGE INDICATOR EXAMPLE ###
n = 10000
my_prperc = mc.progressbar.ProgPercent(n) # 1) init. with number of iterations
for i in range(n):
    # do some computation
    my_prperc.update() 
  Hide output
[100 %] elapsed[sec]: 3.790
Total time elapsed: 3.791 sec

Handling filenames

By default, the folder for generating figures or data is mc.figpath:

In [3]:
print mc.figpath
  Hide output
results/

To generate figures, we assign file names, such as:

In [4]:
name = 'color'
import os
filename = os.path.join(mc.figpath, name)
  Hide output

It is then possible to check if that figures exist:

In [5]:
print filename, mc.check_if_anim_exist(filename), mc.recompute
  Hide output
results/color True False

Importing MayaVi

In [6]:
print mc.import_mayavi.__doc__
  Hide output
None

Mayavi is difficult to compile on some architecture (think Mac OsX), so we allowed the possibility of an ImportError:

In [7]:
print mc.MAYAVI
mc.import_mayavi()
  Hide output
WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.

Import
Imported Mayavi

To force MayaVi to not be imported, set:

In [8]:
# mc.MAYAVI = 'Do NOT Import'
# mc.import_mayavi()
  Hide output

exporting to various formats

It is possible to export motion clouds to many different formats. Here are some examples:

In [9]:
name = 'export'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
z = mc.rectif(mc.random_cloud(mc.envelope_gabor(fx, fy, ft)))
mc.PROGRESS = False
for vext in mc.SUPPORTED_FORMATS:
    print 'Exporting to format: ', vext
    mc.anim_save(z, os.path.join(mc.figpath, name), display=False, vext=vext, verbose=False)
  Hide output
Exporting to format:  .h5
Exporting to format: 
[100 %] elapsed[sec]: 0.298
 .mpg
Saving sequence results/export.mpg
Title: 
                      Started: 10/14/2014 17:41:11
                      Finished: 10/14/2014 17:41:11
                      Total time elapsed: 0.299 sec
                      CPU %: 101.700000
                      Memory %: 2.920246

Total time elapsed: 0.299 sec
[100 %] elapsed[sec]: 0.315

Exporting to format:  .mp4
Saving sequence results/export.mp4
Title: 
                      Started: 10/14/2014 17:41:11
                      Finished: 10/14/2014 17:41:11
                      Total time elapsed: 0.315 sec
                      CPU %: 97.900000
                      Memory %: 2.922750

Total time elapsed: 0.315 sec
[100 %] elapsed[sec]: 0.260

Exporting to format:  .gif
Saving sequence results/export.gif
Title: 
                      Started: 10/14/2014 17:41:11
                      Finished: 10/14/2014 17:41:12
                      Total time elapsed: 0.261 sec
                      CPU %: 102.000000
                      Memory %: 2.923656
Exporting to format: 

Total time elapsed: 0.261 sec
[100 %] elapsed[sec]: 0.272
 .webm
Saving sequence results/export.webm
Title: 
                      Started: 10/14/2014 17:41:12
                      Finished: 10/14/2014 17:41:12
                      Total time elapsed: 0.272 sec
                      CPU %: 102.000000
                      Memory %: 2.923656

Total time elapsed: 0.272 sec
[100 %] elapsed[sec]: 0.291

Exporting to format:  .zip
Saving sequence results/export.zip
Title: 
                      Started: 10/14/2014 17:41:12
                      Finished: 10/14/2014 17:41:13
                      Total time elapsed: 0.292 sec
                      CPU %: 102.200000
                      Memory %: 2.923656
Exporting to format:  .mat
Exporting to format: 

Total time elapsed: 0.292 sec
[100 %] elapsed[sec]: 0.296
 .mkv
Saving sequence results/export.mkv
Title: 
                      Started: 10/14/2014 17:41:15
                      Finished: 10/14/2014 17:41:15
                      Total time elapsed: 0.296 sec
                      CPU %: 102.200000
                      Memory %: 2.934384


Total time elapsed: 0.296 sec

showing a video

To show a video in a notebook, issue:

In [19]:
mc.notebook = True # True by default
mc.in_show_video('export')
  Hide output

Rectifying the contrast

The mc.rectif function allows to rectify the amplitude of luminance values within the whole generated texture between $0$ and $1$:

In [11]:
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
color = mc.envelope_color(fx, fy, ft)
envelope = color * mc.envelope_gabor(fx, fy, ft)
image = mc.random_cloud(envelope)
print 'Min :', image.min(), ', mean: ', image.mean(), ', max: ', image.max()
  Hide output
Min : -581.212000276 , mean:  -8.0823883827e-18 , max:  593.148283647

In [12]:
image = mc.rectif(image)
print 'Min :', image.min(), ', mean: ', image.mean(), ', max: ', image.max()
  Hide output
Min : 0.0590556234669 , mean:  0.5 , max:  0.95

Generating figures

It is easy to generate figures, either directly from a set of parameters for the Motion Clouds (with mc.figure_MC) or from an arbitrary envelope (with mc.figure):

In [17]:
name = 'test_figures_MC'
#initialize
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
mc.figures_MC(fx, fy, ft, name)
mc.in_show_video(name)
  Hide output
/usr/local/lib/python2.7/site-packages/traits/has_traits.py:1541: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  setattr( self, name, value )
[100 %] elapsed[sec]: 0.303
Saving sequence results/test_figures_MC.webm
Title: 
                      Started: 10/14/2014 17:42:17
                      Finished: 10/14/2014 17:42:18
                      Total time elapsed: 0.303 sec
                      CPU %: 102.100000
                      Memory %: 3.831553
test_figures_MC


Total time elapsed: 0.303 sec

In [18]:
name = 'test_figures'
#initialize
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
color = mc.envelope_color(fx, fy, ft)
z = color * mc.envelope_gabor(fx, fy, ft)
mc.figures(z, name)
mc.in_show_video(name)
  Hide output
[100 %] elapsed[sec]: 0.383
Saving sequence results/test_figures.webm
Title: 
                      Started: 10/14/2014 17:42:22
                      Finished: 10/14/2014 17:42:23
                      Total time elapsed: 0.384 sec
                      CPU %: 101.700000
                      Memory %: 5.887246
test_figures


Total time elapsed: 0.384 sec

There are also two specific functions to generate the mayavi vizualizations: cube and visualize.

Testing the different components of the enveloppe

Motion Clouds: testing components of the envelope

In [2]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output
In [3]:
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
#mc.recompute = True
mc.notebook = True
  Hide output

Testing the speed

Here the link to the test page for the component Speed

In [5]:
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
color = mc.envelope_color(fx, fy, ft)
name = 'speed'
z = color*mc.envelope_speed(fx, fy, ft)
mc.figures(z, name)
mc.show_video2(name)
  Hide output
speed

Exploring the orientation component of the envelope around a grating.

Here the link to the test page for the component Grating

In [6]:
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
color = mc.envelope_color(fx, fy, ft)
name = 'grating'
z = color * mc.envelope_gabor(fx, fy, ft)
mc.figures(z, name)
mc.show_video2(name)
  Hide output
grating

Here the link to the test page for the component Radial

In [8]:
B_theta = np.pi/8.
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
name = 'radial'
mc.figures_MC(fx, fy, ft, name, B_theta=B_theta)
verbose = False
mc.show_video2(name)
  Hide output
radial

Testing the color

Here the link to the test page for the component Color

In [10]:
name = 'color'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
z = mc.envelope_color(fx, fy, ft)
mc.figures(z, name)
mc.show_video2(name)
  Hide output
color

Testing competing Motion Clouds with psychopy

Analysis of a Discrimination Experiment

In the psychopy_competition.py script, we implemented an experiment to test whether one could discriminate the dominating motion in a sum of two motion clouds in opposite directions.

Herein, we analyse the data that was collected over different sessions and try to draw some conclusions.

In [27]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

Analysis, version 1

In a first version of the experiment, we only stored the results in a numpy file:

In [28]:
!ls ../data/competing_v1_*npy
  Hide output
../data/competing_v1_bruno_Dec_14_1210.npy  ../data/competing_v1_lup_Dec_12_1013.npy  ../data/competing_v1_meduz_Dec_14_1204.npy
../data/competing_v1_lup_Dec_12_1003.npy    ../data/competing_v1_lup_Dec_14_1201.npy

In [29]:
import glob
fig = plt.figure()
ax = fig.add_subplot(111)
for fn in glob.glob('../data/competing_v1_*npy'):
    results = np.load(fn)
    print 'experiment ', fn, ', # trials=', results.shape[1]
    ax.scatter(results[1, :], results[0, :])
_ = ax.axis([0., 1., -1.1, 1.1])
_ = ax.set_xlabel('contrast')
_ = ax.set_ylabel('perceived direction')
  Hide output
experiment  ../data/competing_v1_bruno_Dec_14_1210.npy , # trials= 50
experiment  ../data/competing_v1_lup_Dec_12_1003.npy , # trials= 50
experiment  ../data/competing_v1_lup_Dec_12_1013.npy , # trials= 50
experiment  ../data/competing_v1_lup_Dec_14_1201.npy , # trials= 10
experiment  ../data/competing_v1_meduz_Dec_14_1204.npy , # trials= 50

In [30]:
alpha = .3
fig = plt.figure(figsize=(12,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
data = []
for fn in glob.glob('../data/competing_v1_*npy'):
    results = np.load(fn)
    data_ = np.empty(results.shape)
    # lower stronger, detected lower = CORRECT is 1
    ax1.scatter(results[1, results[1, :]<.5], 
                1.*(results[0, results[1, :]<.5]==-1), 
                c='r', alpha=alpha)
    # copying data: contrast (from .5 to 1), correct 
    data_[0, results[1, :]<.5] = 1. - results[1, results[1, :]<.5]
    data_[1, results[1, :]<.5] = 1.*(results[0, results[1, :]<.5]==-1)
    # upper stronger, detected lower = INCORRECT is 1
    ax1.scatter(results[1, results[1, :]>.5], 
                1.*(results[0, results[1, :]>.5]==-1), 
                c='b', alpha=alpha)
    # lower stronger, detected upper = INCORRECT is 1
    ax2.scatter(results[1, results[1, :]<.5], 
                1.*(results[0, results[1, :]<.5]==1), 
                c='b', alpha=alpha, marker='x')
    # upper stronger, detected upper = CORRECT is 1
    ax2.scatter(results[1, results[1, :]>.5], 
                1.*(results[0, results[1, :]>.5]==1), 
                c='r', alpha=alpha, marker='x')
    # copying data: contrast (from .5 to 1), correct 
    data_[0, results[1, :]>=.5] = results[1, results[1, :]>=.5]
    data_[1, results[1, :]>=.5] = 1.*(results[0, results[1, :]>=.5]==1)
    data.append(data_)

for ax in [ax1, ax2]:
    ax.axis([0., 1., -.1, 1.1])
    _ = ax.set_xlabel('contrast')
_ = ax1.set_ylabel('detected lower')
_ = ax1.legend(['lower', 'upper'], loc='right')
_ = ax2.set_ylabel('detected upper')
_ = ax2.legend(['lower', 'upper'], loc='right')
  Hide output

(note the subplots are equivalent up to a flip)

One could not fit psychometric curves for this 2AFC. See for instance https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3110332/

let's explore another way:

In [31]:
import numpy as np
import pylab
from scipy.optimize import curve_fit
 
def sigmoid(c, c0, k):
    y = 1 / (1 + np.exp(-k*(c-c0)))
    return y
 
  Hide output
In [33]:
for fn in glob.glob('../data/competing_v1_*npy'):
    results = np.load(fn)
    cdata, ydata = results[1, :], .5*results[0, : ]+.5
    pylab.plot(cdata, ydata, 'o')

    popt, pcov = curve_fit(sigmoid, cdata, ydata)
    print popt
    c = np.linspace(0, 1, 50)
    y = sigmoid(c, *popt)
    pylab.plot(c, y)
    
pylab.ylim(-.05, 1.05)
pylab.legend(loc='best')
  Hide output
[ 0.291  6.067]
[  0.345  13.492]
[  0.287  39.2  ]
[   0.637  160.796]
[  0.362  13.633]

Analysis, version 2

In the second version, we also stored some information about the experiment

In [41]:
from psychopy import misc
mean, slope = [], []
for fn in glob.glob('../data/competing_v2_*npy'):
    results = np.load(fn)
    data = misc.fromFile(fn.replace('npy', 'pickle'))
    print data
    cdata, ydata = results[1, :], .5*results[0, : ]+.5
    pylab.plot(cdata, ydata, 'o')

    popt, pcov = curve_fit(sigmoid, cdata, ydata)
    mean.append(popt[0])
    slope.append(popt[1])
    c = np.linspace(0, 1, 50)
    y = sigmoid(c, *popt)
    pylab.plot(c, y)

pylab.text(0.05, 0.8, 'mean : %0.2f , slope : %0.2f ' %(np.mean(mean), np.mean(slope)))
pylab.ylim(-.05, 1.05)
pylab.legend(loc='best')
  Hide output
{'N_X': 128, 'N_Y': 128, 'observer': u'jean', 'timeStr': 'Sep_03_1536', 'nTrials': 50, 'N_frame': 128, 'N_frame_total': 128, 'screen_width': 2560, 'screen_height': 1440}
{'N_X': 128, 'N_Y': 128, 'observer': u'laurent', 'timeStr': 'Sep_17_1522', 'nTrials': 50, 'N_frame': 128, 'N_frame_total': 128, 'screen_width': 2560, 'screen_height': 1440}

Testing discrimination between Motion Clouds with psychopy

Analysis of a discrimination experiment

In the psychopy_discrimination.py script, we implemented an experiment to test whether the overall shape of motion clouds could influence discrimination of speed.

Herein, we analyse the data that was collected over different sessions and try to draw some conclusions.

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

The experiment

In [2]:
#%pycat psychopy_discrimination.py
  Hide output

Analysis - version 1

In a first version of the experiment, we only stored the results in a numpy file.

In [3]:
!ls  ../data/discriminating_*
  Hide output
../data/discriminating_john_Jan_23_1515.npy  ../data/discriminating_lup_Jan_24_0914.npy  ../data/discriminating_lup_Jan_24_0937.npy
../data/discriminating_lup_Jan_23_1401.npy   ../data/discriminating_lup_Jan_24_0919.npy  ../data/discriminating_v2_lup_Feb_07_1409.pickle
../data/discriminating_lup_Jan_23_1735.npy   ../data/discriminating_lup_Jan_24_0931.npy  ../data/discriminating_v2_lup_Feb_07_1434.pickle

In [4]:
!ls  ../data/
  Hide output
competing.pickle		    competing_v1_meduz_Dec_14_1204.npy	       discriminating_john_Jan_23_1515.npy  discriminating_lup_Jan_24_0931.npy
competing_v1_bruno_Dec_14_1210.npy  competing_v2_anonymous_Feb_07_1558.pickle  discriminating_lup_Jan_23_1401.npy   discriminating_lup_Jan_24_0937.npy
competing_v1_lup_Dec_12_1003.npy    competing_v2_anonymous_Sep_16_1746.pickle  discriminating_lup_Jan_23_1735.npy   discriminating_v2_lup_Feb_07_1409.pickle
competing_v1_lup_Dec_12_1013.npy    competing_v2_jean_Sep_03_1536.npy	       discriminating_lup_Jan_24_0914.npy   discriminating_v2_lup_Feb_07_1434.pickle
competing_v1_lup_Dec_14_1201.npy    competing_v2_jean_Sep_03_1536.pickle       discriminating_lup_Jan_24_0919.npy

In [5]:
import glob
for fn in glob.glob('../data/discriminating_*npy'):
    results = np.load(fn)
    print fn, results.shape
    print('analyzing results')
    print '# of trials :', np.abs(results).sum(axis=1)
    print 'average results: ', (results.sum(axis=1)/np.abs(results).sum(axis=1)*.5+.5)
  Hide output
../data/discriminating_john_Jan_23_1515.npy (2, 50)
analyzing results
# of trials : [ 50.     24.508]
average results:  [ 0.48  1.  ]
../data/discriminating_lup_Jan_23_1401.npy (2, 50)
analyzing results
# of trials : [ 50.     28.126]
average results:  [ 0.66  1.  ]
../data/discriminating_lup_Jan_23_1735.npy (3, 50)
analyzing results
# of trials : [  9.  14.  13.]
average results:  [ 1.  1.  1.]
../data/discriminating_lup_Jan_24_0914.npy (3, 50)
analyzing results
# of trials : [ 17.  21.  12.]
average results:  [ 0.647  0.857  1.   ]
../data/discriminating_lup_Jan_24_0919.npy (3, 50)
analyzing results
# of trials : [ 10.  25.  15.]
average results:  [ 0.7    0.32   0.533]
../data/discriminating_lup_Jan_24_0931.npy (7, 50)
analyzing results
# of trials : [  3.   4.   8.   8.   7.  12.   8.]
average results:  [ 0.667  1.     0.625  0.375  1.     0.167  0.125]
../data/discriminating_lup_Jan_24_0937.npy (7, 50)
analyzing results
# of trials : [  7.   5.   6.  12.  10.   2.   8.]
average results:  [ 0.857  0.8    0.833  0.417  0.2    1.     0.375]

In [6]:
%whos
  Hide output
Variable   Type       Data/Info
-------------------------------
fn         str        ../data/discriminating_lup_Jan_24_0937.npy
glob       module     <module 'glob' from '/usr<...>/lib/python2.7/glob.pyc'>
np         module     <module 'numpy' from '/us<...>ages/numpy/__init__.pyc'>
plt        module     <module 'matplotlib.pyplo<...>s/matplotlib/pyplot.pyc'>
pylab      module     <module 'pylab' from '/us<...>site-packages/pylab.pyc'>
results    ndarray    7x50: 350 elems, type `float64`, 2800 bytes

Another solution is to use:

In [7]:
#data= 1
#np.savez('toto', data=data, results=results)
#print np.load('toto.npz')['data']
  Hide output

or directly psychopy.misc:

In [8]:
from psychopy import misc
  Hide output
In [9]:
#info = misc.fromFile('../data/discriminating.pickle')
  Hide output
In [10]:
alphas = [-1., -.5, 0., 0.5, 1., 1.5, 2.]
fileName = '../data/discriminating_' + info['observer'] + '.pickle'
info['alphas'] = alphas
info['results'] = results
#numpy.savez(fileName, results=results, alphas=alphas)
misc.toFile(fileName, info)
  Hide output
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-1bd6d05d02af> in <module>()
      1 alphas = [-1., -.5, 0., 0.5, 1., 1.5, 2.]
----> 2 fileName = '../data/discriminating_' + info['observer'] + '.pickle'
      3 info['alphas'] = alphas
      4 info['results'] = results
      5 #numpy.savez(fileName, results=results, alphas=alphas)

NameError: name 'info' is not defined

Analysis - version 2

In the second version, we stored more data:

In [11]:
import glob
from psychopy import misc
fig = plt.figure()
ax = fig.add_subplot(111)
for fn in glob.glob('../data/discriminating_v2_*pickle'):
    data = misc.fromFile(fn)
    print fn, results.shape
    print('analyzing results')
    alphas = np.array(data['alphas'])
    print ' alphas = ', alphas
    print '# of trials :', np.abs(data['results']).sum(axis=1)
    print 'average results: ', (data['results'].sum(axis=1)/np.abs(data['results']).sum(axis=1)*.5+.5)
    width = (alphas.max()-alphas.min())/len(alphas)
    ax.bar(data['alphas'] - width/2, (data['results'].sum(axis=1)/np.abs(data['results']).sum(axis=1)*.5+.5), width=width, alpha=.3 )
#    for i_trial in xrange(data['results'].shape[1]):
#        ax.scatter(data['alphas'][data['results']
ax.set_xlabel('alpha')
ax.set_ylabel('% correct')
  Hide output
../data/discriminating_v2_lup_Feb_07_1409.pickle (7, 50)
analyzing results
 alphas =  [-1.  -0.5  0.   0.5  1.   1.5  2. ]
# of trials : [ 10.   4.   5.  11.   6.  11.   3.]
average results:  [ 0.7    1.     0.4    0.545  0.833  0.273  1.   ]
../data/discriminating_v2_lup_Feb_07_1434.pickle (7, 50)
analyzing results
 alphas =  [-1.  -0.5  0.   0.5  1.   1.5  2. ]
# of trials : [  4.   3.  10.  12.   8.   8.   5.]
average results:  [ 0.75   0.667  0.8    0.5    0.125  0.25   0.4  ]

Out[11]:
<matplotlib.text.Text at 0x10d17a610>
In [12]:
# TODO : correct for the number of trials / make global reading / make a new draw of alphas
# TODO : test f_0
  Hide output
In []:
  Hide output

Testing basic functions of Motion Clouds

Motion Clouds: raw principles

Motion Clouds are synthesized textures which aim at having similar characteristics as natural images but with controlled parameters. There are many ways to achieve these results and this notebook aims at showing that different procedures from different communities (neurioscience, modelling, computer vision, ...) may produce similar results.

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

Using Fourier ("official Motion Clouds")

In [2]:
import sys
sys.path.append('..')
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
  Hide output

Using mixtures of images

In [15]:
import scipy.misc
lena = scipy.misc.lena() * 1.
lena -= lena.mean()
lena /= lena.std()
print lena.shape
  Hide output
(512, 512)

In [16]:
plt.imshow(lena, cmap=plt.cm.gray)
  Hide output
Out[16]:
<matplotlib.image.AxesImage at 0x19967c3d0>
In [17]:
lena.shape
  Hide output
Out[17]:
(512, 512)
In [19]:
lena[0, :]
  Hide output
Out[19]:
array([ 0.793,  0.793,  0.793,  0.772,  0.793,  0.689,  0.814,  0.772,
        0.856,  0.772,  0.793,  0.751,  0.647,  0.814,  0.751,  0.647,
        0.689,  0.668,  0.772,  0.772,  0.626,  0.668,  0.626,  0.689,
        0.647,  0.689,  0.647,  0.584,  0.668,  0.626,  0.626,  0.668,
        0.626,  0.709,  0.647,  0.751,  0.709,  0.898,  0.751,  0.877,
        0.877,  0.856,  0.877,  1.002,  0.981,  1.065,  1.023,  0.96 ,
        1.002,  1.002,  0.898,  1.065,  0.918,  0.898,  0.793,  0.772,
        0.501,  0.501,  0.647,  0.354,  0.124, -0.105, -0.126, -0.377,
       -0.565, -0.565, -0.628, -0.67 , -0.753, -0.565, -0.461, -0.586,
       -0.44 , -0.502, -0.419, -0.398, -0.398, -0.398, -0.419, -0.294,
       -0.294, -0.335, -0.356, -0.398, -0.419, -0.314, -0.314, -0.314,
       -0.335, -0.356, -0.356, -0.335, -0.294, -0.314, -0.356, -0.314,
       -0.314, -0.294, -0.294, -0.335, -0.419, -0.377, -0.377, -0.335,
       -0.314, -0.273, -0.147, -0.294, -0.231, -0.126, -0.064, -0.105,
       -0.252, -0.043, -0.064, -0.043,  0.02 , -0.043, -0.064, -0.001,
       -0.001,  0.124,  0.104, -0.022, -0.043,  0.062,  0.166,  0.145,
        0.166,  0.104,  0.229,  0.104,  0.083,  0.145,  0.145,  0.062,
        0.229,  0.25 ,  0.166,  0.208,  0.145,  0.145,  0.083,  0.124,
        0.208,  0.041,  0.104,  0.041,  0.187,  0.124,  0.104,  0.145,
        0.083,  0.166,  0.104,  0.083,  0.229,  0.187,  0.25 ,  0.208,
        0.083,  0.208,  0.124,  0.229,  0.187,  0.187,  0.208,  0.229,
        0.166,  0.229,  0.145,  0.25 ,  0.208,  0.166,  0.187,  0.187,
        0.25 ,  0.187,  0.208,  0.208,  0.166,  0.187,  0.145,  0.166,
        0.187,  0.229,  0.271,  0.229,  0.187,  0.145,  0.187,  0.145,
        0.229,  0.25 ,  0.166,  0.229,  0.25 ,  0.271,  0.187,  0.292,
        0.25 ,  0.208,  0.187,  0.229,  0.229,  0.25 ,  0.166,  0.187,
        0.333,  0.166,  0.145,  0.312,  0.208,  0.145,  0.104,  0.166,
        0.166,  0.062,  0.145,  0.166,  0.166,  0.145,  0.104,  0.166,
        0.145,  0.083,  0.229,  0.187,  0.187,  0.229,  0.145,  0.166,
        0.166,  0.166,  0.145,  0.187,  0.124,  0.145,  0.145,  0.187,
        0.25 ,  0.208,  0.166,  0.229,  0.166,  0.229,  0.208,  0.292,
        0.166,  0.166,  0.187,  0.208,  0.104,  0.104,  0.187,  0.292,
        0.292,  0.375,  0.208,  0.166,  0.166,  0.104,  0.083,  0.083,
        0.166,  0.208,  0.104,  0.145,  0.208,  0.187,  0.145,  0.083,
        0.104,  0.062,  0.145,  0.166,  0.041,  0.124,  0.229,  0.208,
        0.145,  0.062,  0.041,  0.145,  0.062,  0.145,  0.208,  0.062,
        0.041,  0.166,  0.145,  0.083,  0.145,  0.062,  0.145,  0.02 ,
        0.062,  0.124,  0.02 ,  0.062,  0.062,  0.166,  0.041,  0.104,
       -0.064, -0.001, -0.001, -0.064, -0.105, -0.043, -0.064, -0.105,
       -0.147, -0.189, -0.252, -0.231, -0.189, -0.398, -0.461, -0.44 ,
       -0.419, -0.252, -0.168, -0.126, -0.043,  0.166,  0.229,  0.312,
        0.333,  0.417,  0.459,  0.542,  0.605,  0.584,  0.605,  0.73 ,
        0.73 ,  0.73 ,  0.814,  0.709,  0.793,  0.772,  0.647,  0.647,
        0.501,  0.542,  0.605,  0.542,  0.626,  0.563,  0.563,  0.647,
        0.563,  0.689,  0.647,  0.584,  0.689,  0.668,  0.626,  0.626,
        0.73 ,  0.709,  0.647,  0.689,  0.626,  0.605,  0.563,  0.689,
        0.605,  0.563,  0.626,  0.584,  0.626,  0.605,  0.584,  0.605,
        0.605,  0.563,  0.563,  0.709,  0.626,  0.626,  0.584,  0.584,
        0.605,  0.689,  0.626,  0.709,  0.709,  0.689,  0.626,  0.626,
        0.668,  0.73 ,  0.73 ,  0.73 ,  0.689,  0.668,  0.668,  0.584,
        0.647,  0.647,  0.709,  1.107,  1.42 ,  1.545,  1.713,  1.754,
        1.859,  1.838,  1.859,  1.921,  1.963,  1.963,  1.963,  1.963,
        2.005,  1.921,  1.817,  1.671,  1.42 ,  1.002,  0.584, -0.043,
       -0.377, -0.398, -0.482, -0.44 , -0.419, -0.273, -0.398, -0.231,
       -0.168, -0.147, -0.105, -0.147, -0.105, -0.126, -0.064, -0.022,
        0.062, -0.064, -0.105, -0.147, -0.085,  0.02 , -0.064, -0.126,
       -0.022, -0.043,  0.02 ,  0.083, -0.064, -0.105, -0.022, -0.126,
       -0.085, -0.022, -0.022, -0.105, -0.043, -0.064, -0.064, -0.022,
       -0.043, -0.001,  0.104, -0.126, -0.001,  0.041, -0.001,  0.041,
        0.041,  0.062, -0.022, -0.001, -0.001,  0.02 , -0.126, -0.126,
        0.083, -0.043,  0.041,  0.041,  0.041,  0.124, -0.001,  0.041,
        0.02 , -0.064,  0.104, -0.001,  0.041,  0.062, -0.105, -0.022,
       -0.21 , -0.064, -0.168, -0.105, -0.126, -0.189, -0.085, -0.001,
        0.292,  0.835,  0.918,  0.96 ,  0.981,  0.96 ,  0.647,  0.083])
In [5]:
def noise(image=lena):
    for axis in [0, 1]:
        image = np.roll(image, np.random.randint(image.shape[axis]), axis=axis)
    return image
  Hide output
In [6]:
plt.imshow(noise(), cmap=plt.cm.gray)
  Hide output
Out[6]:
<matplotlib.image.AxesImage at 0x117793790>
In [7]:
plt.imshow(noise(), cmap=plt.cm.gray)
  Hide output
Out[7]:
<matplotlib.image.AxesImage at 0x117c95890>

Using ARMA processes

Now, we define the ARMA process as an averaging process with a certain time constant $\tau=30.$ (in frames).

In [8]:
def ARMA(image, tau=30.):
    image = (1 - 1/tau)* image + 1/tau * noise()
    return image
  Hide output

initializing

In [9]:
image = ARMA(lena)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[9]:
<matplotlib.image.AxesImage at 0x118193b50>
In [10]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[10]:
<matplotlib.image.AxesImage at 0x11820bcd0>
In [11]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[11]:
<matplotlib.image.AxesImage at 0x118bccd10>
In [12]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[12]:
<matplotlib.image.AxesImage at 0x11910b8d0>
In [13]:
%cd..
N_frame = 1024
z = np.zeros((lena.shape[0], lena.shape[1], N_frame))
z[:, :, 0] = image
for i_frame in range(1, N_frame): 
    z[:, :, i_frame] = ARMA(z[:, :, i_frame-1])
mc.anim_save(.5 + .5*z, filename='results/arma')
  Hide output
[  0 %] elapsed[sec]: 0.000
/Users/lolo/Dropbox/science/MotionClouds
Saving sequence results/arma.webm
[51150 %] elapsed[sec]: 57.826 | ETA[sec]: -57.713 

Title: 
                      Started: 09/23/2014 11:58:20
                      Finished: 09/23/2014 11:59:18
                      Total time elapsed: 0.000 sec

In [14]:
mc.notebook = True
mc.in_show_video(name='arma')
  Hide output
arma

    Share